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ABSTRACT 

This paper discusses the mathematical 
existence and the numerically-correct 
identification of linear and nonlinear aerodynamic 
impulse response functions. Differences between 
continuous-time and discrete-time system 
theories, which permit the identification and 
efficient use of these functions, will be detailed. 
Important input/output definitions and the 
concept of linear and nonlinear systems with 
memory will also be discussed. It will be shown 
that indicial (step or steady) responses (such as 
Wagner’s function), forced harmonic responses 
(such as Theodorsen’s function or those from 
doublet lattice theory), and responses to random 
inputs (such as gusts) can all be obtained from an 
aerodynamic impulse response function. This 
paper establishes the aerodynamic impulse 
response function as the most fundamental, and, 
therefore, the most computationally efficient, 
aerodynamic function that can be extracted from 
any given discrete-time, aerodynamic system. 
The results presented in this paper help to unify 
the understanding of classical two-dimensional 
continuous-time theories with modem three- 
dimensional, discrete-time theories. First, the 
method is applied to the nonlinear viscous 
Burger’s equation as an example. Next the 
method is applied to a three-dimensional 
aeroelastic model using the CAP-TSD 
(Computational Aeroelasticity Program 
Transonic Small Disturbance) code and then to a 
two-dimensional model using the CFL3D 
Navier-Stokes code. Comparisons of accuracy 
and computational cost savings are presented. 
Because of its mathematical generality, an 
important attribute of this methodology is that it 
is applicable to a wide range of nonlinear, 
discrete-time problems. 


INTRODUCTION 

During the early development of 
mathematical models of unsteady aerodynamic 
responses, the efficiency and power of 
superposition of scaled and shifted fundamental 
responses, or convolution, was quickly 
recognized. This led to the classical Wagner’s 
function 1 , which is the response of a two- 
dimensional airfoil, in incompressible flow, to a 
unit step variation in angle of attack. Similar 
functions such as Kussner’s function, which is 
the response of a two-dimensional airfoil to a 
sharp-edged gust in incompressible flow, were 
developed as well 1 . 

As geometric complexity increased, however, 
the analytical derivation of these time-domain 
fundamental functions became quite complicated 
and, therefore, impractical. Ultimately, 
frequency-domain aerodynamics for three- 
dimensional configurations became the method of 
choice for computing linear unsteady 
aerodynamic responses 2 . For the case where 
geometry- and/or flow-induced nonlinearities are 
significant in the aerodynamic response, time 
integration of the nonlinear equations is 
necessary, as is done in unsteady CFD codes, 
particularly for aeroelastic analyses. As CFD 
codes have grown in complexity and capability, 
there is a very real need to incorporate these codes 
into aeroservoelastic (ASE) analyses, loads 
estimation, and other preliminary design efforts 
in an efficient and accurate manner. Direct 
incorporation of a CFD code into the ASE 
process is currently not practical due to the high 
computational costs and turnaround time required. 
As computational speeds improve and as new 
algorithms are developed to address this problem, 
the practicality of this approach may improve. 
At the moment, however, the efficient 
incorporation of the information provided by a 
CFD code into disciplines such as ASE remains 
a problem. 
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Attempts to address this problem include the 
development of transonic indicial responses 3,4,5 . 
Reference 6 develops models of nonlinear 
aerodynamic maneuvers from an experimental 
database using neural networks. References 7 and 
8 provide reduced-order models for linear and 
linearized solutions about a nonlinear condition. 

In order to develop robust, mathematically- 
correct and efficient nonlinear models of the CFD 
response, a mathematically-formal method is 
required that is well defined in the time and 
frequency domains and that is well defined for 
continuous- and discrete-time systems. The 
discrete-time Volterra theory of nonlinear 
systems fulfills these requirements and was 
applied in the present research. This theory has 
found wide acclaim in the field of nonlinear 
discrete-time systems 9 and nonlinear digital filters 
for telecommunications and image processing 10 , 
to name a small subset of references. 
Applications of this theory to nonlinear, discrete- 
time aerodynamic systems include Tromp and 

Jenkins ', Rodriguez 12 , and Silva 13,14 . In Ref. 13, 
the concept of applying the Volterra theory to the 
development of efficient linear and nonlinear 
aerodynamic impulse responses was presented and 
demonstrated to be feasible for high frequencies. 
In Ref. 14, the identification, and computational 
efficiency of linear discrete-time aerodynamic 
impulse responses, valid for arbitrary inputs, was 
demonstrated using the linear equations within 
the CAP-TSD 15 (Computational Aeroelasticity 
Program - Transonic Small Disturbance) code. 
Nonlinear aerodynamic impulse responses were 
identified using the nonlinear equations within 
the CAP-TSD code but were limited in scope 
because of the particular identification technique 
that was used. The present paper removes these 
limitations by presenting a mathematically- 
correct identification scheme for nonlinear 
responses. Reference 14 represents the first time 
that aerodynamic impulse response functions 
were numerically identified. The concept of 
linear and nonlinear aerodynamic impulse 
response functions introduces a totally new 
perspective on linear and nonlinear, steady and 
unsteady aerodynamics, as will be discussed. 

The purpose of this paper is to introduce 
new, or improved, mathematical developments 
that allow the mathematically-correct 
identification of linear and nonlinear aerodynamic 
impulse responses. The functional classification 
of the discrete-time Navier-Stokes equations that 


enable the correct application of the discrete-time 
Volterra theory to CFD codes is presented. The 
fundamental nature of these responses with 
regards to classical and modem aerodynamic 
theories and the impact of these developments on 
fields such as aeroelasticity and ASE is discussed 
as well. As an illustrative example, the discrete- 
time Volterra theory is applied to the nonlinear 
viscous Burger’s equation. Then the theory is 
applied to a three-dimensional aeroelastic model 
using the CAP-TSD code and then to an airfoil 
in plunge using the CFL3D 16 Navier-Stokes 
code. Comparisons of accuracy and 
computational cost savings are presented. 

MATHEMATICAL PRELIMINARIES 

Discretized Navier-Stokes Equations 

The application of CFD codes involves, in 
general, the application of the discretized Navier- 
Stokes (NS) equations. This is true for the entire 
spectrum of equation levels, from the linear 
equations to the full Navier-Stokes equations, 
including transonic small-disturbance (TSD) and 
Euler equations. The only difference between the 
different equation levels is the number and type 
of simplifying assumptions used to derive the 
resultant governing equations. It is important, 
therefore, to understand the functional nature of 
the NS equations 17 from a mathematical systems 
perspective. 

Upon convergence of an initial, steady-state 
solution, the discretized NS equations form a 
discrete-time, nonlinear, time-invariant system. 
Reynold’s averaging of the NS equations and 
inclusion of turbulence models to provide closure 
does not alter this aspect of the equations. This 
realization, formally stated here for the first time, 
allows the application of techniques routinely 
used in the modeling and design of nonlinear, 
discrete-time filters. In particular, Ref. 18 proves 
that discrete-time, nonlinear, time-invariant 
systems with memory can be modeled arbitrarily 
well using Volterra models, neural networks, or 
radial basis functions. An important attribute of 
Volterra models is that physical interpretation of 
the resulting functions is possible, in the time 
and frequency domains, which often reveals an 
underlying structure of the system. Description 
of the Volterra theory of nonlinear systems is 
presented in Refs. 13 and 14, and the references 
therein. 
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A time-invariant system, also referred to as a 
stationary or autonomous system, is a system 
whose fundamental properties do not change with 
time. An example of a simple, time-invariant, 
nonlinear system is a pendulum. Although the 
full nonlinear equation of a pendulum is certainly 
a function of time which exhibits nonlinear, 
unsteady responses, neither the length of the 
pendulum nor the mass at the end of the 
pendulum are functions of time 19 . 

A time-varying system, also referred to as a 
non-stationary or non-autonomous system, is a 
system whose fundamental properties change 
with time. Fortunately, for many of the 
problems in aircraft unsteady aerodynamics, 
aeroelasticity, and aeroservoelasticity, the 
governing nonlinear equations are time-invariant. 
The linearization of these time-invariant, 
nonlinear equations about an operating point 
yields the familiar time-invariant, linear 
equations that comprise the majority of modern- 
day, linear analyses techniques in these fields. 

The memory of a system, linear or nonlinear, 
is a measure of dependence of the system on 
outputs from previous times. The impulse 
response of a linear system is the “memory” of 
that system. It is a temporal representation of 
the manner in which and the length of time over 
which a unit perturbation remains active in the 
response of the system. Convolution then is 
used to predict the exact response of the linear 
system to an arbitrary input (any and all steady 
and unsteady inputs) because all responses of the 
system are scaled and shifted superpositions of 
this memory function. Likewise, the concept of 
memory functions can be extended to nonlinear 
systems via the Volterra theory of nonlinear 
systems. 

Numerical approximations to ordinary and 
partial differential equations, such as finite- 
difference techniques, are defined by the 
dependence of the response on previous values of 
input and output. Clearly then, time-accurate, 
discretized models, such as finite-difference 
models, are systems with memory, by definition. 
A discretized version of the NS equations (after 
steady- state convergence) is, therefore, a time- 
invariant, nonlinear, discrete-time system and the 
application of the discrete-time Volterra theory to 
this system of equations is a valid mathematical 
approach as proved by Ref. 18. 


Discrete-time Systems 

The modem field of discrete-time signal 
processing 20 is a mathematical systems field that 
addresses substantially more issues than just the 
sampling of a continuous-time signal. A main 
topic in this field is that of digital filter design. 
In digital filter design, there exist mathematical 
concepts that are quite different from their 
continuous-time counterparts. The first of these 
is the unit impulse function, or the Dirac delta 
function. Whereas the continuous-time unit 
impulse is an abstract function, typically 
considered impractical for actual applications 21 or 
sometimes misinterpreted as an indicial (step) 
input 6,22,23 , the discrete-time equivalent, known as 
the unit sample function, is a simple, well- 
defined and extremely useful function. Digital 
filters are designed using this input and its 
resultant output known as the unit sample 
response. The unit sample function is defined as 

u[t] = 1.0 fork=k0 

= 0.0 for k*k0 (1) 


The application of this input to a linear, discrete- 
time system will yield the system’s unit sample 
response, the discrete-time equivalent of the unit 
impulse response. The properties of the unit 
sample response are identical to those of the unit 
impulse response. Both responses completely 
define a linear system and, through convolution, 
the response of the system to any arbitrary input 
can be predicted exactly without actually 
processing the arbitrary input through the 
system. This is because the unit sample 
response captures the system’s complete 
frequency content. 

A linear system’s frequency characteristics 
can be determined by applying multiple sinusoids 
of varying frequency, applying band-limited 
white noise, or by computing the fast Fourier 
transform (FFT) of the unit sample response. 
The application of multiple sinusoids is, 
basically, how linear, frequency-domain, unsteady 
aerodynamics are generated. The band-limited 
white noise technique implies exploration of 
different segments of the system’s bandwidth in a 
piecewise, overlapping, and inefficient fashion. 
The most efficient approach is to compute the 
FFT of the unit sample response, yielding the 
system’s frequency response. This efficiency is 
the result of the fundamental properties of the 
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unit sample response. Additional evidence of 
this efficiency is the fact that the response of the 
system to the multiple sinusoidal inputs and the 
band-limited white noise can be computed via 
convolution of these inputs with the unit sample 
response. Therefore, from the single 
computation of the unit sample response, all 
system responses, from steady (step) to random, 
can be generated as well. This concept is well 
understood and routinely applied in the design of 
digital filters yet appears to be rare in fields 
dominated by continous-time concepts. 

The concept of convolution is another idea 
that is routinely used in digital filter design but 
that is perceived as somewhat abstract, and 
therefore avoided, by the continuous-time 
community. Because it was believed that 
practical application of an impulse to an 
aerodynamic system could not be performed, 
discrete-time aerodynamic impulse responses 
were never identified until recently in Ref. 14. 

Convolution, in discrete-time, is defined as 

y[n] = fh[n-k]x[k] (2) 

k-0 

where h[n-k] is the unit sample response and x[k] 
is the arbitrary input. It is important to 
understand that this is not the discrete-time 
version of Duhamel’s integral 24 , which is the 
convolution of a unit step response with the 
derivative of an arbitrary input. The unit step 
(indicial) response is not the same as the unit 
sample (impulse) response, as some references 
have indicated 6,22,23 . 

The response of a linear system to an 
arbitrary function of time, x[k], can be computed 
via three methods. The first, or trivial method, 
is to process the input through the system itself. 
If the system is complex and computationally 
intensive, significant computational costs, 
including turnaround time, will be incurred. The 
second method is to identify the system’s unit 
step response and then, via convolution with the 
derivative of the arbitrary input, obtain the 
response of the system using 


y[n] = x[0]S[n] + 2S[n-k]x’[k]At (3) 

k=0 

where S[k] is the unit step response and x’[k] is 
the derivative of the arbitrary input. Equation (3) 


is the discrete-time equivalent of Duhamel’s 
integral. The first term in Eq.(3) must, of 
course, be included whenever x[0] is nonzero. 
Equation (3) is the correct discrete-time 
implementation for indicial (or step) 
aerodynamics. It is mathematically-valid if and 
only if the step response is correctly identified 
and applied in Equation (3). The application of 
step functions has typically been a problem in 
computational unsteady aerodynamics because of 
the downwash equation and the perceived problem 
with the derivative of a step input. This issue is 
addressed in a subsequent section of this paper. 

The third method is to identify the system’s 
unit sample response and, via convolution with 
the arbitrary input, x[k], (Eq. (2)), obtain the 
response of the system. Again, proper 
identification of the unit sample response is a 
requirement for the succesful application of this 
method. 

Clearly, for complex and computationally- 
intensive linear systems, the second and third 
methods provide the most efficient method for 
computing responses because repeated execution 
of the system is not required. The unit sample 
response and the unit step response contain all 
the necessary information regarding the system’s 
behavior in a compact form. In addition, the 
derivative of the unit step response is the unit 
sample response so that only one response, the 
step or the unit sample response, is needed to 
compute the other. The derivative of Wagner’s 
function, for example, yields the incompressible, 
aerodynamic impulse response due to plunge for 
a two-dimensional airfoil 25 . Figure 1 was 
obtained using W.P. Jones’ approximation to 
Wagner’s function 21 . Details regarding this result 
and its relation to Theodorsen’s function can be 
found in Ref. 25. 

In this research, the identification and use of 
linear and nonlinear aerodynamic unit sample 
responses is favored over that of the unit step 
responses for the following reasons: ( 1 ) The unit 
step response can be computed via convolution 
of the unit sample response with a step input, 
yielding the steady-state solution; and (2) 
Convolution using the unit sample response 
involves the actual input whereas convolution 
using the unit step response involves the 
derivative of the input, requiring additional, 
unnecessary computational effort. The unit 
sample response is the most compact 

representation of a linear system from which all 
other steady and unsteady responses can be 
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generated. Extension of this concept to nonlinear 
systems then enables the efficient computation of 
nonlinear steady and unsteady responses due to 
arbitrary inputs. 

Identification of linear aerodynamic unit 
sample responses 14 has interesting implications. 
First, it provides an alternative to the forced 
harmonic method for computing unsteady 
aerodynamic forces by computing the unit 
sample responses for each mode and then 
performing the convolutions with sinusoidal 
inputs of varying frequency. This could be done 
more directly by performing a Fourier transform 
of each of the modal unit sample responses. 

The unsteady aerodynamic frequency domain 
may be avoided altogether by performing the 
aeroelastic analyses directly in the time domain 13 . 
This is done by coupling the aerodynamic unit 
sample responses with the linear structure in a 
closed-loop sense and obtaining the time-accurate 
aeroelastic transients. Since the aerodynamic 
unit sample response is valid in the complex 
plane, there is no need for rational function 
approximations 23 (RFAs) that extend the forced- 
harmonic responses, valid only along the 
imaginary axis, to the complex plane via analytic 
continuation. Current methods for generating 
RFAs, limited by a specified frequency range of 
interest to generate a low-order model, are 
actually modeling that portion of the unit sample 
response that contains the particular frequency 
range of interest 13 . The aerodynamic unit sample 
response can also be used to realize a linear, 
discrete-time, state-space system 26 . This 
approach was investigated by the author and will 
be the subject of another paper. 

Linear frequency domain and RFA methods 
are not applicable to nonlinear aerodynamics and, 
consequently, the generation of time- accurate, 
aeroelastic transients is necessary. The discrete- 
time Volterra theory of nonlinear systems, along 
with new mathematical developments presented 
in this paper, provides a formal method for the 
identification of nonlinear unit sample responses. 
This results in significant computational 
efficiency when applied, for example, to a CFD 
code. 

Aerodynamic System Input Definition 

An important conceptual development of Ref. 
14, and its subsequent improvement in the 
present research, was the mathematically-correct 
definition of the input to an unsteady 
aerodynamic system for the discrete-time domain. 


The input function consists of the downwash 
function, which, for the excitation of a given 
mode is written as 

w(x,y,t) = phi’(x,y)m(t) + phi(x,y)*u’(t) (4) 

where phi(x,y) is the modeshape, phi’(x,y) are 
the slopes of the modeshape, u(t) is the 
generalized coordinate, and u’(t) is the derivative 
of the generalized coordinate. The discussion 
will be limited, temporarily, to the linear case. 

The current method for the excitation of 
aeroelastic modes within a CFD code involves 
the definition of a “smooth” function defined as 


u(t) = d 0 *exp(-w(t-t 0 )**2) (5) 

where do is the maximum amplitude desired, w is 
the width, and Xq is the time at which the 
maximum amplitude is reached. This Gaussian 
curve (Equation (5)) is referred to as the 
exponential pulse function. This exponential 
pulse is input to each of the modes of the system 
to obtain the set of exponential pulse responses, 
about a nonlinear steady state solution 2728 that 
are then transformed to the frequency domain for 
use in standard linear analyses techniques. This 
should not be confused with the unit pulse 
response mentioned throughout this paper. 
Whereas the unit pulse input (Eq. (1)) excites all 
the frequencies for a given mode, the exponential 
pulse input will excite only the particular range 
of frequencies defined by the width of the 
exponential pulse. This can be explained using 
Eq. (4) as follows. 

From Equation (5), the downwash equation 
consists of the first term which multiplies u(t) 
by the slopes of the modeshape added to u’(t) 
multiplied by the modeshape. When the shape of 
u(t) is narrowed, then the derivative term, u’(t), is 
much bigger and changes more rapidly than it 
does for the wider pulse, thereby exciting higher 
frequencies. Shape optimization may, therefore, 
have to be performed to obtain the desired 
frequency range of interest. Typically, a ‘"wide” 
pulse is recommended, forcing the u'(t) term to 
be small. 

A critical drawback, however, is that the 
exponential pulse is perceived, erronously, as a 
single input. That is, the fast Fourier Transform 
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(FFT) of the output generalized force is divided 
by the FFT of the perceived single input, u(t), to 
obtain the linearized frequency response function 
for that generalized force. But inspection of Eq. 
(4) clearly shows that the downwash function is a 
two-input function. The user defines u(t) but the 
quantity that is input to the flow solver is Eq. 
(4), which includes the effect of u’(t) as well. 
Because this derivative is computed analytically 
internal to the code, it is invisible to the user, 
although it is certainly not invisible to the flow 
solver. 

Inspection of Eq. (4) for a plunge mode 
reveals that the first term is identically zero 
because the slopes of a plunge mode are zero. 
Therefore, the only temporal function that is 
actually input to the flow solver is u’(t). For a 
plunge mode, the denominator of the frequency 
response function should be the FFT of u’(t), not 
the FFT of u(t). This will be demonstrated using 
convolution with examples from CAP-TSD and 
CFL3D in the results section of this paper. 

It is because this second term of the 
downwash input has been ignored that wiggles 
appear at lower or higher frequencies, depending 
on the input u(t), in early applications of the 
technique 27 . The reason for the success of the 
technique to date is that for most modes, a very 
wide u(t) term results in a very small u’(t) term, 
thereby exciting, predominantly, the lower 
frequency range which is, typically, where most 
analyses are desired anyway. If an accurate 
determination of the entire frequency range of a 
mode is desired, then the second term of the 
downwash function must be included in the 
analysis. In terms of computational efficiency, 
the exponential pulse response does not possess 
any of the mathematical properties of the unit 
sample response nor can it be formally extended 
to nonlinear systems. 

The misinterpretation of the downwash as a 
single input has led to the false conclusion that 
impulse (or unit pulse) and step inputs cannot be 
applied to a CFD code because these inputs will 
result in numerical difficulties. The reasoning is 
that the application of a unit pulse, or unit step, 
input as u(t) would lead to a very large, if not 
infinite, derivative term, u’(t). So typically, a 
step input is modified, or made “smoother”, so 
that the u’(t) does not cause numerical problems. 
These “smoother” responses, however, are not 
mathematically consistent with the strict 
definition of unit pulses or unit step inputs and 
so will yield inaccuracies when used in 


convolution. The unit pulse and unit step 
functions have a very precise mathematical 
description which allows for convolution to be 
applied. Any deviation from this precise 
definition will reduce, or possibly eliminate, the 
accuracy of the convolution. 

Mathematically, the downwash equation (for 
a given mode) is clearly a two-channel input. 
For the linear case, each term of the downwash 
equation can, and should, be treated as a separate 
input channel 14 . For the nonlinear case, the 
response due to the sum of the terms of the 
downwash will not be equal to the sum of the 
separate responses due to each term of the 
downwash. The inputs, however, still need to be 
treated as independent inputs. This difficulty was 
solved by computing a combined unit sample 
response that consists of a unit sample input 
applied to each of the two inputs simultaneously 
while using a deconvolution 25 technique to 
maintain mathematical accuracy. This 
deconvolution technique identified the proper 
temporal function that can be used with the 
combined unit sample response to yield the 
correct final response for the linear case. Since 
the combined motion of the system due to the 
combined inputs of the downwash is the same for 
the linear and nonlinear cases, the same combined 
motion is used in the linear and nonlinear 
convolutions. The effectiveness of this method 
will be presented in the results section of this 
paper. 

VQLTERRA THEORY 

The discrete-time Vol terra series for a 
truncated, second-order, time-invariant, system 
has the form 

y[n] = h 0 + fhj[ n -k]x[k] + 
k=0 

N N 

2 2 h 2tn - kl ,n - k2] x[kl] x[k2] 
kl=0 k.2=0 

( 6 ) 

where y[n] is the response of the nonlinear 
system to x[k], an arbitrary input; ho is the mean 
value about which the response is defined; hi is 
the first-order kernel or the linear unit sample 
response; and h2 is the second-order kernel. 
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Details of the theoretical definitions of this 
method, including identification of the kernels, 
can be found in Refs. 13, 14 and all the 
references therein. As in Refs. 13 and 14, 
modeling of the nonlinear aerodynamic system 
will be limited to identification of the first- and 
second-order kernels. 

An intuitive explanation of the application of 
this approach to a nonlinear system can be stated 
as follows. It is a well-established procedure to 
linearize a nonlinear system by expanding the 
nonlinear terms in a Taylor Series about a chosen 
point. The resultant Taylor Series, if expanded 
to sufficient terms, is an excellent approximation 
to the actual nonlinearity. That is, there are no 
restrictions on its range of applicability regarding 
input amplitudes. As the series is truncated by 
gradual elimination of the higher-order terms 
from highest to lowest, limitations on the range 
of applicability of the series approximation 
become more restrictive until the only term left 
is the linear term, the most severely restricted 
term of all. If higher-order terms are gradually 
added back in to the series approximation, one at 
a time, the accuracy of the approximation is 
improved and the range of applicability is 
increased as well. The present method is, 
therefore, a method that re- in states higher-order 
terms that were removed during the linearization 
of the equations. This will yield improved 
accuracy over the purely linear solution and will 
increase the range of applicability as well. 

Also, when a “small’' (or “linear”) input is 
applied to a nonlinear system, there is an 
implicit assumption of the equivalence between 
the nonlinearity and its series expansion. This is 
evident because it is in the presence of a series 
expansion formulation that a “small” input will, 
in fact, yield the “linear” portion of the response 
since the higher-order terms (second-order and 
above) are much smaller and, therefore, 
negligible. The accepted practice of using a 
“small” amplitude exponential pulse response 
(within a CFD code, for example) to excite only 
the “linear” portion of the response about a 
nonlinear solution implies a series 
approximation to the nonlinearity. As a result, 
this “small” input approach offers additional 
validation to the present application of the 
discrete-time Volterra theory, which seeks to 
identify the next higher-order term after the linear 
term. 

Furthermore, the first-order term is more 
accurate than the purely linear term because the 


first-order term is derived with knowledge of the 
second-order, or higher-order, terms. Therefore, 
for a second-order nonlinearity, the first-order 
term is the proper and correct linearization. The 
first-order term can be considered to represent a 
“mean” value of the response with the second- 
order term representing a higher-order variation 
about that mean. 

The success of linearized aerodynamic 
predictions for certain flight regimes, and under 
certain small perturbation assumptions, is due to 
the fact that highly nonlinear phenomena have a 
negligible impact on the net effect of various 
responses at these conditions. It does not mean 
that rotational, viscous, and turbulent effects 
disappear from the flow at these conditions, but 
rather that these effects do not excite higher-order 
effects sufficiently to affect the overall response. 
Increasing the order of this restricted linearized 
approximation to model higher-order effects is, 
therefore, a logical step. 

The computational efficiency of the present 
technique is due to the following features of the 
method. 1). Identification of the first- and 
second-order kernels eliminates the need to re- 
execute the code. 2). The kernels can be coupled 
with a structure in a closed-loop sense “outside” 
of the CFD code, on a workstation, sidestepping 
the current, very expensive method of solving the 
aeroelastic equations of motion within the CFD 
code. 3). The identification of the kernels is 
geometry independent. The kernel of a three- 
dimensional configuration is, topologically, no 
different from the kernel of a two-dimensional 
configuration. The only difference is the initial 
cost of identification that requires the use of the 
CFD code. The complex CFD model, consisting 
of three spatial variables and one temporal 
variable, is mapped onto the unit sample 
response, a compact function of time only. The 
modal approach and the definition of boundary 
conditions within a CFD code make this 
mapping possible. 4). This technique permits a 
unified approach for generation of compact, 
linearized and nonlinear, steady and unsteady 
models from the same, arbitrarily complex CFD 
model (complete configuration, finest grid, most 
detail), including, of course, stability derivatives. 

RESULTS 

Linear CAP-TSD 

The linear equations within the CAP-TSD 
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code were used for comparisons of unit sample 
and step responses. The computational model is 
a rectangular wing with an aspect ratio of two. 
All results presented are for M=0.9. Shown in 
Figure 2 is a comparison of the plunge unit 
sample response and the plunge unit step 
response. Convolution of the unit sample 
response with a unit step also yields the unit step 
response, as shown in Figure 3. Convolution of 
the unit sample response with the input shown 
in Figure 4, u’(t), yields the exact, CAP-TSD- 
generated result, also shown in Figure 4. 
Convolution of the plunge unit sample response 
with u’(t), instead of u(t), yields the correct 
result, consistent with the discussion regarding 
Equation (4) in a previous section. 

These results demonstrate the relationship 
between a unit sample response and a unit step 
response for a linear unsteady aerodynamic 
system and the connect application of these 
functions. Also, it is important to realize that 
the unit sample response, when convoluted with 
a step input results in the steady-state solution, 
as shown. Therefore, unit sample responses can 
be used for predicting the steady and unsteady 
responses of a system. This applies to the 
nonlinear case as well where the savings in 
computational cost and time are of greater 
importance. 

Viscous Burger’s Equation 

The 1-D viscous Burger’s equation is defined 
as 


du du 

— + u — 
dt dx 



( 7 ) 


and is typically used as a simplified model of the 
Navier-Stokes equations for evaluating the 
effectiveness of numerical methods 29 . It is used 
here to demonstrate the effectiveness of the 
discrete-time Volterra technique. The numerical 
solution is implemented via a simple forward-in- 
time, central-in-space (FTCS) method. 

The identification part of the process consists 
of the generation of the first- and second-order 
kernels of a selected grid point due to 
perturbation of the end-point boundary condition. 
Shown in Figure 5 is the first-order kernel of the 
system, revealing a well-behaved and compact 
function. Shown in Figure 6 are the first twenty 
terms of the symmetric second-order kernel. 
These terms indicate a second-order nonlinear 


memory that goes to zero fairly quickly. 

Shown in Figure 7 is a comparison of several 
responses due to step inputs of increasing 
amplitude for the actual numerical solution, the 
convolution of the first-order kernel with each of 
the inputs, and the convolution of first- and 
second-order kernels with each of the inputs. As 
the amplitude is increased, the error between the 
actual (“true”) response and the first-order 
response increases, indicating an increasing effect 
of the nonlinearity as amplitude is increased. 
Addition of the second-order convolution shows a 
significant improvement in accuracy, as seen in 
Figure 7. The crossing over of the convolved 
response for the largest step response could be an 
indication of a convergence limit or the need for 
additional terms of the second-order kernel. The 
improvement in response with the addition of the 
second-order term is, nontheless, evident. Using 
only the first-and second-order kernels, steady- 
state responses of the nonlinear system can be 
computed without re-execution of the actual 
numerical system. It is interesting to note that, 
for a certain range of amplitudes, the first-order 
response may be sufficient, depending on the 
level of accuracy desired. 

Actual and convolved responses, using the 
same first- and second-order kernels, due to 
sinusoidal inputs were generated 25 . Shown in 
Figure 8 is the comparison for one of these 
inputs. Again, the comparisons were excellent 
with the combined first- and second-order 
response showing the best agreement with the 
actual responses. In the case of a purely linear 
system, these responses could be used to generate 
the frequency response function of the system, as 
is done in the doublet lattice technique for linear 
aerodynamic systems. Therefore, whereas the 
unit sample responses are valid in the complex 
domain, the forced harmonic response, which can 
be generated from the unit sample response, is 
valid only on the imaginary axis. The unit 
sample responses (linear) and first- and second- 
order kernels (nonlinear) do not have any such 
limitation. The only limitation of the nonlinear 
functions is that the radius of convergence of the 
series is limited by the norm of the input 13,14 , 
which depends on the system being investigated. 
These functions are therefore more powerful and, 
at the same time, more efficient than any other 
responses that can be obtained from a given 
system. This is because all other system 
responses are the result of a convolution of the 
system’s unit sample response with some 


8 



arbitrary input. 

Shown in Figure 9 is a comparison of the 
actual, first-order, and first- plus second-oider 
responses due to a quasi-random input from a 
uniform probability distribution. Again, the 
comparisons are reasonable for the first-order 
only and excellent for the first- plus second-oider 
responses. This is analogous to the computation 
of the response of a nonlinear system (aircraft) 
due to a random input, such as a gust. Therefore, 
just as in the linear case, the first- and second- 
order kernels can be used to predict the response 
of the nonlinear system to any arbitrary input, 
which is an infinite set of possible inputs. 

Nonlinear CAP-TSD 

The nonlinear TSD equation is solved within 
the CAP-TSD code for a rectangular wing with a 
NACA0012 airfoil section undergoing plunge 
and an aspect ratio of two at a Mach number of 
0.9. Figure 10 is a comparison of nonlinear 
CAP-TSD responses, due to plunging motions 
of different amplitudes, with the convolved 
results of the first-order kernel with the same 
inputs. The linear CAP-TSD result for the first 
amplitude is also shown for comparison. For 
this mode, the first-order kernel seems to be 
sufficient to capture the range of responses. This 
is not surprising given the nature of the TSD 
equation. The cost for ten of these types of 
responses using CAP-TSD directly is 38,000 
CPU secs and 15 hours turnaround time. The 
cost using the first-order convolution for ten of 
these types of responses is 4,150 CPU secs and 
2.04 hours turnaround time. Most of the cost of 
the first-order convolution is the initial 
identification part of the process since each 
convolution itself took only 75 seconds on a 
workstation. As the need for the response of the 
system to arbitrary inputs (motions) increases, 
the cost of the method decreases because once the 
unit sample responses are obtained, the CFD 
code need not be re-executed. 

Figure 1 1 is a comparison of the actual 
nonlinear CAP-TSD solution for the same wing 
undergoing an arbitrary pitching motion and the 
response obtained by the convolution of the 
combined first-order kernel and the appropriate 
input, obtained as described in an earlier section 
of the paper. The comparison is reasonable, but 
for this mode, the second-order terms are needed 25 . 
The computational efficiency has, however, been 
doubled and is now mathematically correct for 
nonlinear responses. The reason for this is that 


instead of computing two responses per mode 
(one for each term of the down wash function, Eq. 
(4)), only one response per mode is needed. 

CFL3D (version 5.0) 

Navier-Stokes results for a dense-grid RAE 
airfoil 16 with the Spalart-Allmaras turbulence 
model undergoing plunge at M=0.75 were 
computed at a time step of 0.001. The RAE 
airfoil grid is presented in Figure 12. 

Comparison of the CFL3D responses with 
the first-order convolved responses, as well as a 
linear response, are shown in Figure 13. The 
comparisons are excellent, with decreasing 
accuracy as the amplitude increases, similar to 
the viscous Burger’s equation results. As 
amplitude increases, so does the need for second- 
order kernels. Details for this case and higher 
Mach numbers (increased nonlinearity) can be 
found in Ref. 25. 

These results prove the applicability of 
discrete-time, nonlinear, unit sample responses at 
the NS equation level, as discussed in the 
beginning of the paper. 

The cost of each CFL3D run was about 
2,000 CPU seconds. The cost of the first-order 
kernel identification was 400 CPU seconds 
because the kernel goes to zero in less than 100 
time steps. The cost of each convolution, 
performed on a workstation, was 30 seconds. 
The most important point, however, is that a 
compact model has been identified that is valid 
for a range of amplitudes without re-execution of 
the code. 

CONCLUSIONS 

The mathematically correct and numerically- 
accurate identification of linear and nonlinear, 
discrete-time aerodynamic impulse responses was 
presented. For the linear case, the aerodynamic 
impulse response functions were used to 
reproduce exactly the responses of a linearized 
three-dimensional aeroelastic CFD model to 
arbitrary aeroelastic input motions at a fraction of 
the computational cost and time. It was shown 
that the response to step (steady), sinusoidal, and 
random inputs can all be computed from an 
impulse response function, establishing the 
aerodynamic impulse response function as the 
most fundamental aerodynamic function that can 
be extracted from a discrete-time, aerodynamic 
system. 
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For the nonlinear case, the existence, 
identification, and application of nonlinear, 
discrete- time, aerodynamic impulse responses was 
presented. Applications of the method to the 
nonlinear viscous Burger’s equation revealed the 
existence of well-behaved first- and second-order 
impulse response functions. The method was 
then applied to nonlinear aeroelastic CFD models 
using the CAP-TSD and CFL3D codes. The 
results prove the existence of these functions for 
complex, three-dimensional CFD models and 
their application demonstrates their accuracy and 
computational efficiency. 
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Semi-chord Lengths of Travel 

Figure 1 Impulse response due to unit plunge, derived from 
W.P. Jones' approximation to Wagner's function. 
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Figure 4 Comparison of responses, due to arbitrary plunging 
motion, linear, M=0.9, DT=0.001 
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Figure 2 Comparison of unit sample and indicial responses 
linear, M=0.9, DT=0.001 
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Figure 5 The first-order kernel of the response in velocity 

to unit perturbation for the viscous Burger's equation. 
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Figure 3 Comparison of indicial responses, linear, 
M=0.9, DT=0.001 
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Burger's equation. 
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Figure 7 Comparison of actual (A), first-order (1), and 

first-+ second-order (1+2) step responses for the 
viscous Burger's equation. 
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Figure 10 Comparison of actual nonlinear and first-order 
convolved for three plunging motions, 

M=0.9, DT=0.001 
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Figure 8 Comparison of actual (A), first-order (1), and 

first— + second-order (1+2) harmonic responses for 
the viscous Burger's equation. 
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Figure 1 1 Comparison of linear and nonlinear, actual and 

convolved, responses for CAP-TSD model pitching 
at M=0.9. 
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Figure 9 Comparison of actual (A), first-order (1), and 

first-+ second-order (1+2) responses due to quasi- 
random imput for the viscous Burger's equation. 
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Figure 12 RAE 2822 airfoil grid, from Ref. 16. 
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